stitches in R
IMPORTANT NOTES
- Necessary to avoid R variable names like
variable.namesince.in python is somewhat analogous to$in R. - Follow the instructions in the
stitches-in-R-setupR markdown so that this notebook will work. If you’ve followed that markdown,stitchesshould be installed in ther-reticulatevirtual environment for this notebook and should be callable.
Setup
R
#### R BLOCK ####
library(reticulate)## Warning: package 'reticulate' was built under R version 3.5.2# # indicate that we want to use a specific virtualenv
use_virtualenv("r-reticulate", required =TRUE)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(ggplot2)
# time slice for the map plotting:
map_time <- '2100-12-31'python
#### Python Block ####
import stitches
import pandas as pd
import pkg_resources
import xarray as xr
import numpy as np
pd.set_option('display.max_columns', None)
# The CMIP6 ESM we want to emulate and the variables we want to
# emulate
esm = ['CanESM5']
vars1 = ['tas', 'pr', 'hurs', 'psl']Load Hector data
pulled a quick rcp4.5 from HectorUI
#### R BLOCK ####
hector_tgav <- read.csv('data/Hector-data-2022-09-22.csv', skip=3, stringsAsFactors = FALSE)
hector_exp <- unique(hector_tgav$scenario)
head(hector_tgav)## scenario year variable value units
## 1 RCP Standard-RCP-4.5 1800 Tgav 0.1176974 degC
## 2 RCP Standard-RCP-4.5 1801 Tgav 0.1212491 degC
## 3 RCP Standard-RCP-4.5 1802 Tgav 0.1243894 degC
## 4 RCP Standard-RCP-4.5 1803 Tgav 0.1272345 degC
## 5 RCP Standard-RCP-4.5 1804 Tgav 0.1298589 degC
## 6 RCP Standard-RCP-4.5 1805 Tgav 0.1323124 degCReshape the Hector data to be stitches target
First we reshape the Hector time series to match the CMIP6 ESM GSAT smooth anomaly time series that STITCHES operates on: NOTE we should probably add this as a function to stitches NOTE Hector time series start in 1800 - there’s no science reason we couldn’t do matching over that window and have a gridded netcdf of an extra 50 years of history.
#### R BLOCK ####
hector_tgav %>%
mutate(variable = 'tas',
ensemble = 'hectorUI1',# doesn't affect the matching or stitching
model = 'Hector') %>%
rename(experiment = scenario) %>%
filter(year >= 1850) %>%
select(variable, experiment, ensemble, model, year, value, -units) ->
target_tgav
print(head(target_tgav))## variable experiment ensemble model year value
## 1 tas RCP Standard-RCP-4.5 hectorUI1 Hector 1850 0.08854785
## 2 tas RCP Standard-RCP-4.5 hectorUI1 Hector 1851 0.09338570
## 3 tas RCP Standard-RCP-4.5 hectorUI1 Hector 1852 0.09959573
## 4 tas RCP Standard-RCP-4.5 hectorUI1 Hector 1853 0.10672673
## 5 tas RCP Standard-RCP-4.5 hectorUI1 Hector 1854 0.11386913
## 6 tas RCP Standard-RCP-4.5 hectorUI1 Hector 1855 0.10890346Convert Hector Tgav to stitches windows
Now that we have shaped it as stitches is expecting, we can use stitches functions to calculate the matching windows for this target data
#### Python Block ####
target_data = stitches.fx_processing.get_chunk_info(
stitches.fx_processing.chunk_ts(df = r.target_tgav, n=9))## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:121: FutureWarning: Passing a set as an indexer is deprecated and will raise in a future version. Use a list instead.
## extra_info = df[extra_columns].drop_duplicates()
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:163: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
## fx_dx_info = fx_dx_info.append(row)
## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_processing.py:171: SettingWithCopyWarning:
## A value is trying to be set on a copy of a slice from a DataFrame.
## Try using .loc[row_indexer,col_indexer] = value instead
##
## See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
## extra_info["i"] = 1print(target_data.head())## experiment model ensemble variable start_yr end_yr year \
## 0 RCP Standard-RCP-4.5 Hector hectorUI1 tas 1850 1858 1854
## 1 RCP Standard-RCP-4.5 Hector hectorUI1 tas 1859 1867 1863
## 2 RCP Standard-RCP-4.5 Hector hectorUI1 tas 1868 1876 1872
## 3 RCP Standard-RCP-4.5 Hector hectorUI1 tas 1877 1885 1881
## 4 RCP Standard-RCP-4.5 Hector hectorUI1 tas 1886 1894 1890
##
## fx dx
## 0 0.113869 -0.011713
## 1 0.067891 0.011509
## 2 0.141379 0.002686
## 3 0.162569 -0.029602
## 4 -0.045458 0.014252Load the stitches archive data for matching
## ensemble variable model experiment start_yr end_yr year \
## 6996 r10i1p1f1 tas CanESM5 ssp126 1850 1858 1854
## 6997 r10i1p1f1 tas CanESM5 ssp126 1859 1867 1863
## 6998 r10i1p1f1 tas CanESM5 ssp126 1868 1876 1872
## 6999 r10i1p1f1 tas CanESM5 ssp126 1877 1885 1881
## 7000 r10i1p1f1 tas CanESM5 ssp126 1886 1894 1890
##
## fx dx
## 6996 -1.360399 0.016472
## 6997 -1.297474 -0.001508
## 6998 -1.200856 0.010901
## 6999 -1.218503 -0.016847
## 7000 -1.285527 0.006000Match
#### python block ####
my_recipes = stitches.make_recipe(target_data, archive_data,
non_tas_variables=['pr', 'psl', 'hurs'],
tol = 0.13, N_matches = 10000, reproducible = True)## You have requested more recipes than possible for at least one target trajectories, returning what canprint(my_recipes.head())## target_start_yr target_end_yr archive_experiment archive_variable \
## 0 1850 1858 historical tas
## 1 1859 1867 historical tas
## 2 1868 1876 historical tas
## 3 1877 1885 historical tas
## 4 1886 1894 historical tas
##
## archive_model archive_ensemble stitching_id \
## 0 CanESM5 r24i1p1f1 RCP Standard-RCP-4.5~hectorUI1~1
## 1 CanESM5 r14i1p1f1 RCP Standard-RCP-4.5~hectorUI1~1
## 2 CanESM5 r14i1p1f1 RCP Standard-RCP-4.5~hectorUI1~1
## 3 CanESM5 r24i1p1f1 RCP Standard-RCP-4.5~hectorUI1~1
## 4 CanESM5 r1i1p1f1 RCP Standard-RCP-4.5~hectorUI1~1
##
## archive_start_yr archive_end_yr \
## 0 2003 2011
## 1 2003 2011
## 2 2003 2011
## 3 2003 2011
## 4 2003 2011
##
## tas_file \
## 0 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 1 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 2 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 3 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 4 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
##
## hurs_file \
## 0 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 1 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 2 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 3 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 4 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
##
## pr_file \
## 0 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 1 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 2 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 3 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 4 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
##
## psl_file
## 0 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 1 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 2 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 3 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...
## 4 gs://cmip6/CMIP6/CMIP/CCCma/CanESM5/historical...Unless otherwise specified by the argument
reproducible = True, the draws stitches takes from the pool
of matches is stochastic.
## You have requested more recipes than possible for at least one target trajectories, returning what canWithin each of these ensembles, there’s no envelope collapse. If we were to concatenate them together into a super-ensemble, there would be. So it’s not as powerful as a very large ensemble of collapse free realizations but it’s still a distinct ensemble to consider.
stitch Tgav
#### python block ####
stitched_global_temp = stitches.gmat_stitching(my_recipes)## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_stitch.py:344: FutureWarning: In a future version of pandas, a length 1 tuple will be returned when iterating over a groupby with a grouper equal to a list of length 1. Don't supply a list with a single grouper to avoid this warning.
## for name, match in rp.groupby(['stitching_id']):print(stitched_global_temp.head())## year value variable stitching_id
## 0 1850 0.106926 tas RCP Standard-RCP-4.5~hectorUI1~1
## 1 1851 0.154738 tas RCP Standard-RCP-4.5~hectorUI1~1
## 2 1852 0.118152 tas RCP Standard-RCP-4.5~hectorUI1~1
## 3 1853 0.094706 tas RCP Standard-RCP-4.5~hectorUI1~1
## 4 1854 0.024321 tas RCP Standard-RCP-4.5~hectorUI1~1## /Users/snyd535/.virtualenvs/r-reticulate/lib/python3.9/site-packages/stitches/fx_stitch.py:344: FutureWarning: In a future version of pandas, a length 1 tuple will be returned when iterating over a groupby with a grouper equal to a list of length 1. Don't supply a list with a single grouper to avoid this warning.
## for name, match in rp.groupby(['stitching_id']):plot Tgav
if you’re really pressed for computing time, you can pair this with pattern scaling and get the mean spatial field. But we can do internal variability at the gridded multivariate scale, so we do.
stitch gridded
Using the above recipes, stitches can produce new gridded netcdf files for each recipe, for multiple variables. Creating multiple netcdfs of monthly data is slower than focusing on GSAT. Still much faster than an ESM, and once the recipes are set, very parallelizable. But for tractability, we are just going to stitch the gridded netcdfs for one ensemble member.
#### python block ####
# Just do one realization
recipe = my_recipes[my_recipes['stitching_id'] == 'RCP Standard-RCP-4.5~hectorUI1~1'].reset_index(drop=True).copy()We must specify a directory for the netcdfs to be saved to. We will use simply do it directly here. On my computer, constructing these 4 netcdfs took ~2-3 minutes.
#### python block ####
stitches.gridded_stitching(out_dir='data', rp =recipe)## ['Stitching gridded netcdf for: CanESM5 tas RCP Standard-RCP-4.5~hectorUI1~1']
## ['data/stitched_CanESM5_tas_RCP Standard-RCP-4.5~hectorUI1~1.nc', 'data/stitched_CanESM5_hurs_RCP Standard-RCP-4.5~hectorUI1~1.nc', 'data/stitched_CanESM5_pr_RCP Standard-RCP-4.5~hectorUI1~1.nc', 'data/stitched_CanESM5_psl_RCP Standard-RCP-4.5~hectorUI1~1.nc']plot gridded
- we will do maps of all 4 variables at a specified month and year
- we will do a time series of each variable in a specific grid cell
## [1] "data/stitched_CanESM5_hurs_RCP Standard-RCP-4.5~hectorUI1~1.nc"
## [2] "data/stitched_CanESM5_pr_RCP Standard-RCP-4.5~hectorUI1~1.nc"
## [3] "data/stitched_CanESM5_psl_RCP Standard-RCP-4.5~hectorUI1~1.nc"
## [4] "data/stitched_CanESM5_tas_RCP Standard-RCP-4.5~hectorUI1~1.nc"python maps and time series
It’s a lot easier to use python to load in the netcdfs and select a grid cell or a time slice, save those off as simpler data frames than a full netcdf, and plot in R. The reshaping set of code that we often do in R to work with netcdf data is just so much faster if we’ve pre-sliced the data to a single grid cell or single time step in python.
tas
TODO just make it a function
#### python block ####
py_files = pd.DataFrame(r.files, columns = ['file_name'])
# load the netcdf:
f = py_files.loc[py_files['file_name'].str.contains('tas')].values[0]
tas_nc = xr.open_dataset(f[0])
# time slice map:
tas_time_slice = tas_nc.sel(time = r.map_time).copy()
# time series for single grid cell:
# lon and lat values for a grid cell near College Park, MD, home of JGCRI:
cp_lat = 38.9897
cp_lon = 180 + 76.9378 # this is probably not actually right
# lat and lon coordinates closest
abslat = np.abs(tas_nc.lat - cp_lat)
abslon = np.abs(tas_nc.lon-cp_lon)
c = np.maximum(abslon, abslat)
([lon_loc], [lat_loc]) = np.where(c == np.min(c))
lon_grid = tas_nc.lon[lon_loc]
lat_grid = tas_nc.lat[lat_loc]
tas_grid_slice = tas_nc.sel(lon = lon_grid, lat=lat_grid,
time = slice('1850-01-01', '2099-12-31')).copy()
del(tas_nc)Take a look at this data in R:
#### R block ####
# map data:
print(py$tas_time_slice$tas)## <xarray.DataArray 'tas' (lat: 64, lon: 128)>
## array([[256.15045, 255.97285, 255.76503, ..., 256.75378, 256.55612, 256.3434 ],
## [255.52583, 255.14264, 254.7989 , ..., 256.9645 , 256.42502, 255.96776],
## [254.48996, 254.06548, 253.68169, ..., 256.22208, 255.54071, 254.98186],
## ...,
## [272.51633, 273.26175, 274.00467, ..., 269.78442, 270.852 , 271.72083],
## [272.34247, 272.59592, 272.86353, ..., 271.3565 , 271.73337, 272.0625 ],
## [272.1695 , 272.23404, 272.2921 , ..., 271.9601 , 272.0346 , 272.099 ]],
## dtype=float32)
## Coordinates:
## time datetime64[ns] 2100-12-31
## * lat (lat) float64 -87.86 -85.1 -82.31 -79.53 ... 79.53 82.31 85.1 87.86
## * lon (lon) float64 0.0 2.812 5.625 8.438 ... 348.8 351.6 354.4 357.2
## Attributes:
## units: K
## variable: tas
## experiment: historical
## ensemble: r24i1p1f1
## model: CanESM5
## stitching_id: RCP Standard-RCP-4.5~hectorUI1~1
## recipe_location: data/stitched_CanESM5_tas_RCP Standard-RCP-4.5~hectorU...tas_val <- py$tas_time_slice$tas$values
tas_lon <- py$tas_time_slice$lon$values
tas_lat <- py$tas_time_slice$lat$values
grid <- expand.grid(list(lon=tas_lon,lat=tas_lat))
tas <- cbind(grid,
value=t(matrix(aperm(tas_val,c(2,1)), 1,length(tas_lat)*length(tas_lon))))
as.character(py$tas_time_slice$time) %>%
substr(., 37, 43) ->
tas_time
ggplot(tas, aes(x = lon, y = lat, color=value, fill = value)) + geom_tile() +
ggtitle(paste(tas_time, 'tas (K)'))
# grid time series data:
print(py$tas_grid_slice$tas)## <xarray.DataArray 'tas' (time: 3000)>
## array([267.57404, 267.1224 , 273.28598, ..., 289.8179 , 277.07382, 273.7582 ],
## dtype=float32)
## Coordinates:
## * time (time) datetime64[ns] 1850-01-31 1850-02-28 ... 2099-12-31
## lat float64 37.67
## lon float64 255.9
## Attributes:
## units: K
## variable: tas
## experiment: historical
## ensemble: r24i1p1f1
## model: CanESM5
## stitching_id: RCP Standard-RCP-4.5~hectorUI1~1
## recipe_location: data/stitched_CanESM5_tas_RCP Standard-RCP-4.5~hectorU...cp_tas <- data.frame(value = py$tas_grid_slice$tas$values)
data.frame(time = as.character(py$tas_grid_slice$time$values)) %>%
separate(time, into = c('year', 'month', 'day'), sep = '-') %>%
select(-day) %>%
mutate(year = as.integer(year),
month = as.integer(month),
time_index = as.integer(row.names(.))) ->
cp_time
ggplot(cbind(cp_time, cp_tas), aes(x = time_index, y = value)) + geom_line() +
ggtitle('College Park monthly tas, 1850-2100 (K)') 
ggplot(cbind(cp_time, cp_tas) %>% filter(year >= 2021), aes(x = time_index, y = value)) + geom_line() +
ggtitle('College Park monthly tas, 2021-2100 (K)') 
pr
py_files = pd.DataFrame(r.files, columns = ['file_name'])
# load the netcdf:
f = py_files.loc[py_files['file_name'].str.contains('pr')].values[0]
pr_nc = xr.open_dataset(f[0])
# time slice map:
pr_time_slice = pr_nc.sel(time = r.map_time).copy()
# time series for single grid cell:
# lon and lat values for a grid cell near College Park, MD, home of JGCRI:
cp_lat = 38.9897
cp_lon = 180 + 76.9378 # this is probably not actually right
# lat and lon coordinates closest
abslat = np.abs(pr_nc.lat - cp_lat)
abslon = np.abs(pr_nc.lon-cp_lon)
c = np.maximum(abslon, abslat)
([lon_loc], [lat_loc]) = np.where(c == np.min(c))
lon_grid = pr_nc.lon[lon_loc]
lat_grid = pr_nc.lat[lat_loc]
pr_grid_slice = pr_nc.sel(lon = lon_grid, lat=lat_grid,
time = slice('1850-01-01', '2099-12-31')).copy()
del(pr_nc)Take a look at this data in R:
# map data:
print(py$pr_time_slice$pr)## <xarray.DataArray 'pr' (lat: 64, lon: 128)>
## array([[3.245999e-06, 3.262219e-06, 3.159829e-06, ..., 3.286227e-06,
## 3.319206e-06, 3.363079e-06],
## [2.650309e-06, 2.566143e-06, 2.707141e-06, ..., 2.846019e-06,
## 2.849012e-06, 2.591466e-06],
## [8.970674e-07, 7.867402e-07, 6.428324e-07, ..., 1.386523e-06,
## 1.146150e-06, 1.092138e-06],
## ...,
## [1.823617e-05, 1.563593e-05, 1.459936e-05, ..., 2.877988e-05,
## 2.629005e-05, 2.359155e-05],
## [1.739398e-05, 1.668454e-05, 1.647657e-05, ..., 2.011145e-05,
## 1.950161e-05, 1.870858e-05],
## [1.419007e-05, 1.645207e-05, 1.535799e-05, ..., 1.454292e-05,
## 1.429739e-05, 1.409253e-05]], dtype=float32)
## Coordinates:
## time datetime64[ns] 2100-12-31
## * lat (lat) float64 -87.86 -85.1 -82.31 -79.53 ... 79.53 82.31 85.1 87.86
## * lon (lon) float64 0.0 2.812 5.625 8.438 ... 348.8 351.6 354.4 357.2
## Attributes:
## units: kg m-2 s-1
## variable: pr
## experiment: historical
## ensemble: r24i1p1f1
## model: CanESM5
## stitching_id: RCP Standard-RCP-4.5~hectorUI1~1
## recipe_location: data/stitched_CanESM5_pr_RCP Standard-RCP-4.5~hectorUI...pr_val <- py$pr_time_slice$pr$values
pr_lon <- py$pr_time_slice$lon$values
pr_lat <- py$pr_time_slice$lat$values
grid <- expand.grid(list(lon=pr_lon,lat=pr_lat))
pr <- cbind(grid,
value=t(matrix(aperm(pr_val,c(2,1)), 1,length(pr_lat)*length(pr_lon))))
as.character(py$pr_time_slice$time) %>%
substr(., 37, 43) ->
pr_time
ggplot(pr, aes(x = lon, y = lat, color=value, fill = value)) + geom_tile() +
ggtitle(paste(pr_time, 'pr (kg m-2 s-1)'))
# grid time series data:
print(py$pr_grid_slice$pr)## <xarray.DataArray 'pr' (time: 3000)>
## array([8.414188e-06, 8.546974e-06, 5.744913e-06, ..., 6.263057e-06,
## 1.439180e-05, 2.652666e-06], dtype=float32)
## Coordinates:
## * time (time) datetime64[ns] 1850-01-31 1850-02-28 ... 2099-12-31
## lat float64 37.67
## lon float64 255.9
## Attributes:
## units: kg m-2 s-1
## variable: pr
## experiment: historical
## ensemble: r24i1p1f1
## model: CanESM5
## stitching_id: RCP Standard-RCP-4.5~hectorUI1~1
## recipe_location: data/stitched_CanESM5_pr_RCP Standard-RCP-4.5~hectorUI...cp_pr <- data.frame(value = py$pr_grid_slice$pr$values)
data.frame(time = as.character(py$pr_grid_slice$time$values)) %>%
separate(time, into = c('year', 'month', 'day'), sep = '-') %>%
select(-day) %>%
mutate(year = as.integer(year),
month = as.integer(month),
time_index = as.integer(row.names(.))) ->
cp_time
ggplot(cbind(cp_time, cp_pr), aes(x = time_index, y = value)) + geom_line() +
ggtitle('College Park monthly pr (kg m-2 s-1), 1850-2100') 
ggplot(cbind(cp_time, cp_pr) %>% filter(year >= 2021), aes(x = time_index, y = value)) + geom_line() +
ggtitle('College Park monthly pr (kg m-2 s-1), 2021-2100') 
hurs
py_files = pd.DataFrame(r.files, columns = ['file_name'])
# load the netcdf:
f = py_files.loc[py_files['file_name'].str.contains('hurs')].values[0]
hurs_nc = xr.open_dataset(f[0])
# time slice map:
hurs_time_slice = hurs_nc.sel(time = r.map_time).copy()
# time series for single grid cell:
# lon and lat values for a grid cell near College Park, MD, home of JGCRI:
cp_lat = 38.9897
cp_lon = 180 + 76.9378 # this is probably not actually right
# lat and lon coordinates closest
abslat = np.abs(hurs_nc.lat - cp_lat)
abslon = np.abs(hurs_nc.lon-cp_lon)
c = np.maximum(abslon, abslat)
([lon_loc], [lat_loc]) = np.where(c == np.min(c))
lon_grid = hurs_nc.lon[lon_loc]
lat_grid = hurs_nc.lat[lat_loc]
hurs_grid_slice = hurs_nc.sel(lon = lon_grid, lat=lat_grid,
time = slice('1850-01-01', '2099-12-31')).copy()
del(hurs_nc)Take a look at this data in R:
# map data:
print(py$hurs_time_slice$hurs)## <xarray.DataArray 'hurs' (lat: 64, lon: 128)>
## array([[87.71517 , 87.8501 , 87.895195, ..., 87.50938 , 87.51608 , 87.50802 ],
## [87.6163 , 87.74478 , 87.97272 , ..., 86.70137 , 86.976845, 87.31027 ],
## [87.53369 , 88.48645 , 89.40988 , ..., 84.7092 , 85.60706 , 86.5044 ],
## ...,
## [94.53843 , 94.14827 , 94.08032 , ..., 96.15185 , 95.80203 , 95.10472 ],
## [95.57695 , 95.3057 , 95.14706 , ..., 96.73186 , 96.430595, 95.97026 ],
## [95.035416, 95.10253 , 95.101585, ..., 95.09781 , 95.09777 , 95.04169 ]],
## dtype=float32)
## Coordinates:
## time datetime64[ns] 2100-12-31
## * lat (lat) float64 -87.86 -85.1 -82.31 -79.53 ... 79.53 82.31 85.1 87.86
## * lon (lon) float64 0.0 2.812 5.625 8.438 ... 348.8 351.6 354.4 357.2
## Attributes:
## units: %
## variable: hurs
## experiment: historical
## ensemble: r24i1p1f1
## model: CanESM5
## stitching_id: RCP Standard-RCP-4.5~hectorUI1~1
## recipe_location: data/stitched_CanESM5_hurs_RCP Standard-RCP-4.5~hector...hurs_val <- py$hurs_time_slice$hurs$values
hurs_lon <- py$hurs_time_slice$lon$values
hurs_lat <- py$hurs_time_slice$lat$values
grid <- expand.grid(list(lon=pr_lon,lat=pr_lat))
hurs <- cbind(grid,
value=t(matrix(aperm(hurs_val,c(2,1)), 1,length(hurs_lat)*length(hurs_lon))))
as.character(py$hurs_time_slice$time) %>%
substr(., 37, 43) ->
hurs_time
ggplot(hurs, aes(x = lon, y = lat, color=value, fill = value)) + geom_tile() +
ggtitle(paste(hurs_time, 'hurs (%)'))
# grid time series data:
print(py$hurs_grid_slice$hurs)## <xarray.DataArray 'hurs' (time: 3000)>
## array([59.83445 , 61.397877, 58.244843, ..., 39.531998, 70.94262 , 62.537365],
## dtype=float32)
## Coordinates:
## * time (time) datetime64[ns] 1850-01-31 1850-02-28 ... 2099-12-31
## lat float64 37.67
## lon float64 255.9
## Attributes:
## units: %
## variable: hurs
## experiment: historical
## ensemble: r24i1p1f1
## model: CanESM5
## stitching_id: RCP Standard-RCP-4.5~hectorUI1~1
## recipe_location: data/stitched_CanESM5_hurs_RCP Standard-RCP-4.5~hector...cp_hurs <- data.frame(value = py$hurs_grid_slice$hurs$values)
data.frame(time = as.character(py$hurs_grid_slice$time$values)) %>%
separate(time, into = c('year', 'month', 'day'), sep = '-') %>%
select(-day) %>%
mutate(year = as.integer(year),
month = as.integer(month),
time_index = as.integer(row.names(.))) ->
cp_time
ggplot(cbind(cp_time, cp_hurs), aes(x = time_index, y = value)) + geom_line() +
ggtitle('College Park monthly hurs (%), 1850-2100') 
ggplot(cbind(cp_time, cp_hurs) %>% filter(year >= 2021), aes(x = time_index, y = value)) + geom_line() +
ggtitle('College Park monthly hurs (%), 2021-2100') 
psl
sea level pressure can’t have a CP time series
py_files = pd.DataFrame(r.files, columns = ['file_name'])
# load the netcdf:
f = py_files.loc[py_files['file_name'].str.contains('psl')].values[0]
psl_nc = xr.open_dataset(f[0])
# time slice map:
psl_time_slice = psl_nc.sel(time = r.map_time).copy()
del(psl_nc)Take a look at this data in R:
# map data:
print(py$psl_time_slice$psl)## <xarray.DataArray 'psl' (lat: 64, lon: 128)>
## array([[ 99825.836, 99850.695, 99875.39 , ..., 99752.64 , 99776.56 ,
## 99801.05 ],
## [ 99993.89 , 100066.93 , 100133.33 , ..., 99746.07 , 99831.99 ,
## 99915.15 ],
## [100046.87 , 100116.68 , 100172.305, ..., 99757.41 , 99865.12 ,
## 99962.6 ],
## ...,
## [100074.8 , 100086.89 , 100104.336, ..., 100081.53 , 100071.18 ,
## 100069.2 ],
## [100266.875, 100287.08 , 100309.414, ..., 100222.4 , 100234.23 ,
## 100249.15 ],
## [100575.07 , 100585.805, 100597.19 , ..., 100547.25 , 100555.75 ,
## 100565.04 ]], dtype=float32)
## Coordinates:
## time datetime64[ns] 2100-12-31
## * lat (lat) float64 -87.86 -85.1 -82.31 -79.53 ... 79.53 82.31 85.1 87.86
## * lon (lon) float64 0.0 2.812 5.625 8.438 ... 348.8 351.6 354.4 357.2
## Attributes:
## units: Pa
## variable: psl
## experiment: historical
## ensemble: r24i1p1f1
## model: CanESM5
## stitching_id: RCP Standard-RCP-4.5~hectorUI1~1
## recipe_location: data/stitched_CanESM5_psl_RCP Standard-RCP-4.5~hectorU...psl_val <- py$psl_time_slice$psl$values
psl_lon <- py$psl_time_slice$lon$values
psl_lat <- py$psl_time_slice$lat$values
grid <- expand.grid(list(lon=psl_lon,lat=psl_lat))
psl <- cbind(grid,
value=t(matrix(aperm(psl_val,c(2,1)), 1,length(psl_lat)*length(psl_lon))))
as.character(py$psl_time_slice$time) %>%
substr(., 37, 43) ->
pr_time
ggplot(psl, aes(x = lon, y = lat, color=value, fill = value)) + geom_tile() +
ggtitle(paste(pr_time, 'psl (Pa)'))